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Abstract 

The model of laminated wave turbulence puts forth a novel com- 
putational problem - construction of fast algorithms for finding exact 
solutions of Diophantine equations in integers of order 10 12 and more. 
The equations to be solved in integers are resonant conditions for 
nonlinearly interacting waves and their form is defined by the wave 
dispersion. It is established that for the most common dispersion as 
an arbitrary function of a wave-vector length two different generic 
algorithms are necessary: (1) one-class-case algorithm for waves in- 
teracting through scales, and (2) two-class-case algorithm for waves 
interacting through phases. In our previous paper we described the 
one-class-case generic algorithm and in our present paper we present 
the two-class-case generic algorithm. 

PACS numbers: 47.10.-g, 47.27.De, 47.27.T, 02.60.Pn 
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1 INTRODUCTION 



The most general problem setting of the wave turbulence theory can be 
regarded in the form of a nonlinear partial differential equation 



where £ and M denote linear and nonlinear part of the equation correspond- 
ingly, the linear part has the standard wave solutions of the form 



where the amplitude A may depend on space variables but not on time, and 
a small parameter e shows that only resonant wave interactions are taken 
into account. The dispersion function uj = u>(k) can be easily found by 
substitution of ip into the linear part of the initial PDE, C((p) = 0, while 
d t «-> ioJ and d Xs <-> ik s , and resonance conditions have the form 



for n interacting waves with wave- vectors fcj, i = 1, 2, n. For most physical 
applications it is enough to regard n = 3 or n = 4, and the most common 
form of dispersion function is 



(for instance, capillary, gravitational and surface water waves, planetary 
waves in the ocean, drift waves in tokamak plasma, etc.) 

The model of laminated wave turbulencepQ describes two co-existing lay- 
ers of turbulence - continuous and discrete - which are presented by real 
and integer solutions of Sys.Q correspondingly. The continuous layer is 
described by classical statistical methods [2] while for the discrete layer new 
algorithms have to be developed. It was shown in [3] that an arbitrary integer 
lattice (m, n), each node of the lattice denoting a wave- vector k = (m,n), 
can be divided into some clusters (classes) and there are two types of solu- 
tions of Sys.Q: those belonging to the same class and those belonging to 
different classes. Mathematically, a class is described as a set of wave- vectors 
for which the values of the dispersion function have the same irrationality. 
For instance, if the dispersion function has the form u = \/m 2 + n 2 , then a 
class is described as follows: 



ip = A exp % [kx — ut] 
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where 7 is a natural number and q is a square-free integer. Physically, it 
means that waves are interacting over the scales, that is, each two interact- 
ing waves generate a wave with a wavelength different from the wave lengths 
of the two initial waves. Interactions between the waves of different classes 
do not generate new wavelengths but new phases. 

In our preceding paper j3] we presented a generic algorithm for computing 
all integer solutions of Sys. (Op) within one class. Four- wave interactions among 
2-dimensional gravitational water waves were taken as the main example, in 
this case Sys. (CD) takes form: 

(ml + nf) 1 ^ + (ml + n^) 1 ^ = (m\ + nl) 1 ^ + (m\ + n 2 ^) 1 ^ 

m l + m 2 = m 3 + m A (2) 

n% + n 2 = n 3 + n 4 

and classes are defined as Cl q = {^q}, where q, called class index, are 
all natural numbers containing every prime factor in degree smaller 4 and 7, 
called weight, all natural numbers. It can be proven that if all 4 wave- vectors 
constructing a solution of Sys. ([2]) do not belong to the same class, then the 
only possible situation is following: all the vectors belong to two different 
classes Cl qi , Cl q2 and the first equation of Sys.© can be rewritten then as 

71^1+72^2 = 71^1+72^2 (3) 

with some 71,72 € N and q±, qi being class indexes. In the present paper we 
deal with this two-class case. 



2 COMPUTATIONAL PRELIMINARIES 

As in the previous paper pE], we are going to find all solutions of Eq.(j2]) in 
some finite domain D, i.e. \rrii\, \rii\ < D for some D G N. The first case 
has been studied for D = 1000, where 7r c /(10 3 ) = 384145 classes have been 
encountered. The straightforward approach, not making use of classes, con- 
sumes, as for the first case, at least 0(D 5 ) operations and is out of question 
(see jl], Sec. 3. 2.1 for discussion of this point). 

Straightforward application of classes also does not bring much. The 
Eq.([n]) is now trivial - but classes are interlocked through linear conditions. 
Even if for each pair of classes we could detect interlocking and find solu- 
tions, if any, in 0(1) operations (which is probably the case, though we did 
not prove it), the overall computational complexity is at least 7r c i(D) 2 - i.e. 
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not much less than 0(D A ). For D = 1000 this implies 1.5 ■ 10 11 pairwise class 
matches which is outside any reasonable computational complexity limits. 

The trouble with this approach - as, for that matter, with virtually any 
algorithm consuming much more computation time than the volume of its 
input and output data implies - is, that we perform a lot of intermediary 
calculations, later discarded. We develop an algorithm performing every 
calculation with a given item of input data just once (or a constant number 
of times). First of all we notice that Eq.([3]) can be rewritten as 

(m 2 1L + n 2 1L ) 1/4 = {m\ R + n\ R ) 1 

< Kl + 4) 1/4 = ( m iR + n i R ) 1 

miL - m 1R = -m 2L + m 2R 
niL - Uxr = -u 2 l + n 2R 

V 

where qi,q 2 are two different class indexes and 71,72 - the corresponding 
weights. 

Definition. For any two decompositions of a number ^yfq into a sum of two 
squares (see fllj)) the value 5 m = mi J — m R is called m- deficiency, S n = riL — n R 
is called n-deficiency and 8 m>n = (S m , 5 n ) - deficiency point. 

We immediately see that for two interacting waves their deficiencies must 
be equal: 5 lm = m 1L - m 1R = -m 2L + m 2R = 5 2m , 5 ln = n 1L - n 1R = 
—n 2 L + n 2 R = 5 2n . For a given weight 7, every two decompositions of 7 4 g into 
a sum of two squares yield, in general, four deficiency points with S m , 5 n > 0. 
Consider unsigned decompositions mi,mfi,ni,% >> 0. Assuming > 
n^R, < nR the four points are (mi + uir, n L + n R ), (mi + iur, —n L + 
Ur), {m L — m R , n L — n R ), {m L — ttir, —n L + n R ) and four (symmetrical) 
points in each of the other three quadrants of the (m, n) plane. 

Definition. The set of all deficiency points of a class for a given weight, 
A^, is called its ^-deficiency set. The set of all deficiency points of a class, 
A q , is called its deficiency set. 

The objects defined above play the main role in our algorithm, so we 
compute as an illustrative example the for the number 50. The number 50 
has three decompositions into sum of two squares, namely, 50 = l 2 + 7 2 = 
5 2 + 5 2 = 7 2 + l 2 , and nonnegative deficiency points of decomposition pairs 



71^1 
l2^q~2 



(4) 
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are (5, 5; 7,1), (5, 5; 1,7), (1,7; 7,1). They constitute a subset of the 
deficiency set Ag , namely, the 12 points with m > 0,n > 0. In each of 
three other quadrants of the (m, n, ) plane there lie 36 more points of this 
set, symmetrical to the ones shown with respect to the coordinate axes. 

The crucial idea behind the algorithm of this paper is very simple and 
follows immediately from the exposition above: 

Sys.(TJ]) has a solution with vectors belonging to the two different 
classes Cl qi ,Cl q2 if and only if their deficiency sets have a non-void 
intersection, A 9i nA 92 ^ 0, i.e. some elements belong to both classes. 



3 ALGORITHM DESCRIPTION 

Calculation of relevant class indexes q by a sieve-like procedure, admissible 
weights 7 and decomposition of 7 4 g into sum of two squares have all been 
treated in full detail in [I]. One new feature we introduced here is, that 
immediately after generating the array of class bases q we purge away those 
which, whatever the admissible weight 7, do not have a decomposition into 
a a sum of two squares 7 4 g = m 2 + n 2 with both m < D, n < D. For 
the problem considered in [I] this would be superfluous because virtually all 
these classes have been anyhow filtered away according to another criterium 
(M(q) = 1, Dec(q) < 4) which does not apply here. In this way we exclude 
100562 classes from the 384145 which the sieving procedure returns. 

Evidently for any deficiency point <5 mi „ inequalities \S m \ < 2D, \S n \ < 
2D hold. And if deficiency sets of two classes have a non-void intersection, 
they also have an intersection over points with non- negative \8 m \, \8 n \. So we 
start with declaring a two-dimensional array arDeficiency(0..2D,0..2D) of 
type byte which serves storing and processing deficiency sets of the classes. 
The array is initialized with all zeroes. 

3.1 The Five-Pass Procedure 
3.1.1 Pass 1: Marking deficiency points 

In the first pass for every class q in the main domain D we generate its 
deficiency set D q . Notice that after generating deficiency set of the class for 
each weight 7 and uniting them we must check for doubles and eventually 
get rid of them. Next, for every deficiency point (5 m , 5 n ) of the class we 
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increment the value of the corresponding element of the array by 1, except 
elements with value 255 whose values are not changed. 

3.1.2 Pass 2: Discarding non-interacting classes 

In the second pass we generate deficiency sets once more and for every point 
of the deficiency set of a class check the values of the corresponding point 
of ar Deficiency. If all these values are equal to 1, no waves of the class 
participate in resonant interactions and the class is discarded from further 
considerations. 

For the problem considered this pass excludes just a few (313) classes, 
so the time gain is very modest. However, we include this step into the 
presentation for two reasons. First of all, it had to be done as no possibility 
of reducing the number of classes considered as much as possible and as 
soon as possible (before the most time-consuming steps) may be neglected. 
Second, though giving not much gain for solution of the problem at hand, 
this elimination techniques may play a major role in further applications of 
our algorithm. 

3.1.3 Pass 3: Linking interaction points to interacting vectors 

In the third pass we generate a more detailed deficiency set for each class, 
i.e. for all classes not discarded in the previous pass: for every deficiency 

— * 

point 5 rritn we store q, 7, mi, ul, tur, ur. We do not discard duplicates as 
we did in the previous two passes. Then we revisit the corresponding points 
of ar Deficiency and to each point whose value is larger than 1 link the 
structure (q, 7, thl, Ul, m,R, ur) described above. 

3.1.4 Pass 4: Gathering interaction points 

In the fourth pass we go through the array ar Deficiency once more and store 
every point with value greater than one in an array ar De f 'iciency S 'ol(l . .2D , 0..1). 
We also relink structures linked to deficiency points to corresponding points 
of the new array. 

3.1.5 Pass 5: Extracting solutions 

— * 

The four passes above leave us with an array of points 5 m ^ n and to each of 
these points a list of structures (q l , 7*, m l L , n l L , m R , n R ) is linked (no less than 
two different q l ). In general, a linked list is here most appropriate. Every 
combination of two structures linked to the same point and having different 
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q % yields a solution of Sys.©. From every solution found, we obtain four 
solutions changing signs of m l ,n l in the general case, i.e. for m l ,n l nonzero. 

Notice that theoretically we could skip Pass 4 and extract solutions di- 
rectly from the array arDeficiency. However, this is not reasonable for 
implementation reasons, and Pass 4 is not very time-consuming. 

3.1.6 Implementation remarks. 

Implementing the algorithm described above, we took a few language-specific 
shortcuts that will be briefly described here. 

Passes 1 and 2 have been implemented one-to-one as described above. 
However, manipulating linked lists in VBA involves considerable overhead 
and for the problem considered in this paper we do not need the complete 
functionality of linked lists, i.e. inserting into/deleting from intermediate 
positions of the list. Our main data structure for Pass 3-5 is a simple two- 
dimensional array arSolHalves(l..A r AfAr ( i e j, 0..7) and in a single line of this 
array we store: 

• the class base q; 

• the coordinates of deficiency point d m ,d n ; 

• the coordinates of two wave vectors belonging to this deficiency point; 

• the index in the array arSolHalves of the next line belonging to the 
same deficiency point. 

which is demonstrated in FigJH below. Here, NMNdef is the number of m,n- 




Figure 1: Array simulation of deficiency point lists: the overall ar- 
ray/list structure 
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vectors linked to all deficiency points to which vectors belonging to two or 
more classes belong (6692832 for D = 1000). We generate the deficiency set 
of each class and fill all members of this line of the array except the last one 
in the process, deficiency point by deficiency point. The last member is filled 
later and in the following way. 

For this pass we also declare an auxiliary array arDeficiencesPrev(1..2D, 
1..2D) initialized with zeroes. Having added a new line {q, d m , d n , ttiil, nn, mm, nm, 
to arSolHalves, we look up the value indd m ,d n of arDeficiencesPrev(<i m , d n ). If 
it is zero (this deficiency point being visited the first time) we just assign this 
point the value of the index of the new line in the array arSolHalves. Oth- 
erwise we first assign arSolHalves(mG?d m ,<z n > 7) the value of the current line's 
index in arSolHalves, then write this number to arDeficiencesPrev(rf m , d n ) 
(see FigJ2D- 



Array of previous data 
fields' indexes 



Index in array of solution 
halves 



Data fields in array 
of solution halves 




Figure 2: Array simulation of deficiency point lists: data fields in 
detail. 

A numerical example for this procedure is given in Table 1. In this way, 
the array index of the next "list" member is stored with the previous one, 
except evidently the last one, where the corresponding field stays zero. 
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Index 


Q 






mi 




uir 


n R 


Nextlndex 


1 


1 


1 


1 





1 


1 





117 


117 


1 


1 


1 


-119 


120 


120 


-119 


1241 


1241 


4 


1 


1 


-1 


2 


2 


-1 


2921 


2921 


8 


1 


1 


-2 


3 


3 


-2 


4958 


4958 


12 


1 


1 


-3 


4 


4 


-3 


8107 


8107 


19 


1 


1 


-4 


5 


5 


-4 


10304 




















6692782 


273559 


1 


1 


-995 


996 


996 


-995 


6692802 


6692802 


273567 


1 


1 


-996 


997 


997 


-996 


6692816 


6692816 


273575 


1 


1 


-997 


998 


998 


-997 


6692828 


6692828 


273580 


1 


1 


-998 


999 


999 


-998 


6692832 


6692832 


273583 


1 


1 


-999 


1000 


1000 


-999 






Table 1. A few lines of the table containing solution halves for the deficiency 
point d = (1, 1) (beginning and end of the sequence). 



3.2 Computational Complexity 

Consider computational complexity of these steps. 

3.2.1 Pass 1 

For a single class index q and weight 7, generating deficiency points in the 
first step consumes < 0(log 2 (7 4 g)) operations because every number X has 
no more than 0(logX) decompositions into two squares which we combine 
pairwise to find deficiency points. Decompositions themselves can be found 
in 0(log(7 4 g)) time [5]. There are (_D/g) 1//4 admissible weights to class in- 
dex q, so the overall complexity for a class can be estimated from above as 
log 2 DD 1//4 . Merging deficiency points into A q can be done in O(XlogX) 
time for number of points X, i.e. no more than 0(log 2 DD 1 ^ log(log 2 DD 1 ^)) = 
Oilog'DD 1 ^) 

Taking a rough upper estimate for the number of classes 0(D 2 ), we obtain 
an estimate 0(log 3 DD 2 25 ). Incrementing the points of arDeficiency is 
linear on the point number of the set A q and need not be considered for 
computational complexity separately. 
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3.2.2 Pass 2 



The same complexity estimate holds for the second pass. Notice that, hav- 
ing enough memory, or using partial data loading/unloading similar to that 
used in [6] , we could preserve deficiency sets calculated on the first pass and 
not recalculate them here. However, this would not significally improve the 
overall computational complexity of the algorithm. 

We can not give an a priori estimate for the number of classes discarded 
at the second pass, so we ignore it and hold the initial rough upper estimate 
0(D 2 ) for the number of classes in our further considerations. 

3.2.3 Pass 3 

In the third pass, to every point 5 m>n (no more than 0(log 2 DD 1 ^) of them) 
we link the values (q\ Y, m ii n L> m% Ri n% Ft) f° r "which this point has been struck. 
This, as well as linking to the points of ar Deficiency is, clearly, linear on 
the number of points and does not raise the computational complexity. 

3.2.4 Pass 4 

Complexity of the fourth pass can be estimated as follows. Suppose the 
worst case, i.e. no classes are discarded at step 2 and every deficiency point 
is a solution point, i.e. for every 5 m , n = (S m , 5 n ) no less than two classes 
have deficiency points with the same d m , d n . Then we must make no more 
than 0(log 2 DD 225 ) entries into the new array ar Deficiency Sol. We must 
relink no more than the mean of 0(log 2 D) structures per point, which gives 
an upper estimate of 0(log 4 DD 225 ) time for the pass. However, remember 
that the estimate for the deficiency point number has been made on the 
assumption that all (m l L , n l L , m l R , n l R ) generate distinct deficiency points. In 
simple words, for every point linked to X > 1 structures we obtain X — 1 
less solution points. Now elementary consideration allow us to improve the 
estimate to 0(log 2 DD 2 25 ) time. 

3.2.5 Pass 5 

We did not manage to obtain a reasonable estimate for the computational 
complexity of the fifth step. For the worst case of all structures grouped at 
a single point, the estimate is 0(log 4 DD 45 ) - but this is not realistic. If 
the number of solution points is 0(log 2 DD 2 25 ) and the number of linked 
deficiencies is bounded by some number c, then we can make an estimate 
0(c 2 log 2 DD 2 - 25 ). This, however, is also not quite the case as our numerical 
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simulations show). However, this last step deals with solution extraction and 
extracts them in linear time per solution. Any algorithm solving the problem 
has to extract solutions, so we can be sure that our step 5 is optimal - even 
without any estimate of its computational complexity. Summing up, we ob- 
tain the overall upper estimate of computational complexity 0(log 3 DD 225 ) 
reached at steps 1 and 2 plus the time needed for solution extraction. 

4 DISCUSSION 

Our algorithm has been implemented in the VBA programming language; 
computation time (without disk output of solutions found) on a low-end PC 
(800 MHz Pentium III, 512 MB RAM) is about 10 minutes. Some overall 
numerical data is given in the two Figures below. The number of solutions 
for the 2-class-case depending on the partial domain is shown in the Figj3l 
Both curves are almost ideal cubic lines. Very probably they are cubic lines 
asymptotically - the question is presently under study. 
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Figure 3: Number of solutions in partial domains m^rij < D (curve 
marked diamonds) and m 2 + n 2 < D 2 (curve marked circles), D = 
50, 100,. ..1000. 

Partial domains chosen in FigJH]are of two types: squares < D, just 

for simplicity of computations, and circles m 2 + n 2 < D 2 , more reasonable 
choice from physical point of view (in each circle all the wave lengths are 
< D). The curves in the Figj3] are very close to each other in the domain 
D < 500 though number of integer nodes in a corresponding square is D 2 
and in a circle with radius D there are only ttD 2 /4 integer nodes. This indi- 
cates a very interesting physical phenomenon: most part of the solutions is 
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constructed with the wave vectors parallel and close to either axe X or axe Y. 



On the other hand, the number of solutions in rings (L>-50) 2 < mf+n 2 < 
D 2 (corresponds to the wavelengths between D — 50 and D) grows nearly 
perfectly linearly. Of course the number of solutions in a circle is not equal to 
the sum of solutions in its rings: a solution lies in some ring if and only if all 
its four vectors lie in that same ring. That is, studying solutions in the rings 
only, one excludes automatically a lot of solutions containing vectors with 
substantially different wave lengths simultaneously, for example, with wave 
vectors from the rings D — 50 and D + 100. This "cut" sets of solutions can 
be of use for interpreting of the results of laboratory experiments performed 
for investigation of waves with some given wave lengths (or frequencies) only. 




Figure 4: The multiplicities histogram. 

Another important characteristic of the structure of the solution set is 
multiplicity of a vector which describes how many times a given vector is a 
part of some solution. The multiplicity histogram is shown in FigJH On the 
axis X the multiplicity of a vector is shown and on the axis Y the number of 
vectors with a given multiplicity. The histogram of multiplicities is presented 
in Fig. HI it has been cut off - multiplicities go very high, indeed the vector 
(1000,1000) takes part in 11075 solutions. 

Similar histograms computed for different 1-class-cases show that most 
part of the vectors, 70 — 90%% for different types of waves, take part in one 
solution, e.g. they have multiplicity 1. This means that triads or quartets 
are, so to say, the "primary elements" of a wave system and we can explain its 
most important energetic characteristics in terms of these primary elements. 
The number of vectors with larger multiplicities decreases exponentially when 
multiplicity is growing. The very interesting fact in the 2-class-case is the 
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existence of some initial interval of small multiplicities, from 1 to 10, with 
very small number of corresponding vectors. For instance, there exist only 7 
vectors with multiplicity 2. Beginning with multiplicity 11, the histogram is 
similar to that in the 1-class-case. 

This form of the histogram is quite unexpected and demonstrates once 
more the specifics of the 2-class-case compared to the 1-class-case. As one 
can see from the multiplicity diagram in Fig. HI the major part in 2-class- 
case is played by much larger groups of waves with the number of elements 
being of order 40: each solution consists of 4 vectors, groups contain at 
least one vector with multiplicity 11 though some of them can take part in 
the same solution. This sort of primary elements can be a manifestation of 
a very interesting physical phenomenon which should be investigated later: 
triads and quartets as primary elements demonstrate periodic behavior and 
therefore the whole wave system can be regarded as a quasi-periodic one. On 
the other hand, larger groups of waves may have chaotic behavior and, being 
primary elements, define quite different way of energy transfer through the 
whole wave spectrum. 
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